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Abstract. The Hubbard-I approximation is generalized to allow for direct evaluation 
of the equal-time anomalous two-electron propagator for Hubbard model on two- 
dimensional square lattice. This propagator is compared against the quantum Monte 
Carlo data obtained by Aimi and Imada [J. Phys. Soc. Jpn. 76, 113708 (2007)] in the 
limit of strong electron-electron interaction. The Hubbard-I predictions are in a good 
qualitative agreement with the Monte Carlo results. In particular, d-wave correlations 
decay as cr~ 3 ("free electron" behaviour), if separation r exceeds 2-3 lattice constants. 
However, the Hubbard-I approximation underestimates coefficient c by a factor of 
about three. We conclude that the Hubbard-I approximation, despite its simplicity 
and artefacts, captures the qualitative behaviour of the two-particle propagator for 
the Hubbard model, at least for moderate values of r. 
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1. Introduction 

1.1. General outlook 

Consider the Hubbard model's Hamiltonian p]: 

H = H kin + H int - n 2j n ia , (1) 

i,<r 

where 

#kin = ~t ^ C L C i+a<7> #fot = U Yl n ^ n »- ( 2 ) 
i,a,CT i 

Here c| CT and c icr are the creation and annihilation operators for an electron with spin 
projection a on site i of the square lattice, t is the hoping integral, U is the on-site 
Coulomb repulsion, a are vectors connecting the site i with its nearest neighbors, fi 
is chemical potential, n- 1<T = c\ a c ia is the operator of electron number, and a means 
noi-sigma. 

This Hamiltonian is commonly used to describe systems of strongly correlated 
electrons. Among these, we mention such important examples as high-Tc 
superconducting cuprates, manganites, cobaltites, etc. However, there is no commonly 
accepted method that allows to solve Hubbard model for a generic value of doping. 

One of the most intriguing problems of the Hubbard model is the question 
if the model captures the essential physics of the high-temperature copper-oxide 
superconductors. At low interaction one expects that the model exhibits Kohn-Luttinger 
superconductivity [21 El HI [5] . Unfortunately, the critical temperature in such a regime 
is very small. Thus, it is necessary to study the Hubbard model at moderate or 
large interaction strengths. Such a regime is tractable only when the concentration 
of electrons is low [61 [71 [8j. However, the electron concentration for all known copper- 
oxide superconductors is close to one electron per unit cell. Little rigorous knowledge 
is available for this limit. Different approaches return quite controversial answers: 
some Monte Carlo (MC) studies argue against superconductivity [9j, or against high-T c 
superconductivity [10] in Hubbard model, other numerical works favor existence of the 
superconducting ground state with high condensation energy [TT] . 

In the moderate and strong coupling limits there is no guidance from analytical 
studies as well: in this parameter region the majority of theoretical methods become 
unreliable. Thus, one is forced to utilize uncontrollable devices with unknown accuracy. 

In such a situation it appears useful to investigate reliability of uncontrollable 
schemes. A particularly convenient uncontrollable approach for the analysis of the 
Hubbard model is Hubbard-I approximation pQ. Hubbard-I, being invalid for U < zt 
(where z is the number of the nearest neighbor sites on the lattice), could be considered 
as a good first approximation in the case of the strong correlations U 3> zt. Its evident 
advantage is its simplicity. The method allows deriving analytical results in many cases 
or performing numerical calculations at low computational cost. 
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The clarity and transparency of the Hubbard-I approximation explains its 
popularity among the researchers working in the field of strongly correlated electronic 
systems. For us, it is especially important that, after the discovery of high-Tc 
superconductors, the method has been applied to the superconductivity as well. Namely, 
Hubbard-I approximation is used for study of superconductivity in the Hubbard or 
extended Hubbard model [T2], [13], in the Hubbard model with phonons [HI EEl [16], in 
the Hubbard model with attractive interaction [TTJIH]. A related approach, the two-pole 
approximation (which, in some respects, is superior to Hubbard-I approximation), and 
its generalization are applied to study possible superconductivity of the usual Hubbard 
model [IH] and of the multi-band Hubbard model [20]. Similar techniques are used in 
[21| [22J. Systematic expansion in orders of t/U is developed in [23]. It contains the 
Hubbard-I approximation special case. 

1.2. Our results 

Despite broad use of the Hubbard-I approximation, its ability to account for 
superconducting properties is untested. In this paper we address this issue by comparing 
the Cooper pair propagator found with the help of generalized Hubbard-I method 
against the same propagator calculated within the MC framework (in this paper we 
refer to our calculations of the two-particle Green's function as the generalized Hubbard- 
I approximation to distinguish it from the classical Hubbard-I scheme [1] designed for 
evaluation of the single-particle Green's function). 

Such comparison has became possible only recently due to significant progress of 
computational techniques [HHEI]. We use the MC data of Aimi and Imada [ID] , who 
analyze whether the Hubbard model is enough to capture the phenomenon of high-Tc 
superconductivity. These authors use an advanced numerical approach (pre- projected 
Gaussian-basis Monte Carlo). It evaluates the equal-time two-particle propagator in 
d x 2 _ y 2 -channel Pd(r) without any prior assumptions about the structure of the ground 
state wave function. The MC propagator on 10x10 square lattice does not exhibit 
superconducting correlations up to (7 = It. Instead, when r grows, Pdij) decays 
algebraically until finite size effects set in. 

On our part, we generalize the classic Hubbard-I scheme [I] and derive the equation 
of motion for the Cooper pair propagator. We determine Pd for the doped Hubbard 
model ([I]). It is shown that at large U/t the results predicted by the Hubbard-I approach 
are qualitatively consistent with numerical data obtained in Ref. [ID] . In particular, 
both approaches agree that the d x 2_ y 2 -wave pairing correlations decay as cr -3 ("free 
electron" behaviour). The Hubbard-I approximation underestimates coefficient c by a 
factor of order three. Our findings may be viewed as an attempt to establish the limits 
of applicability of the Hubbard-I scheme. 

The paper is organized as follows. In Sect. [2] we derive the equation of motion for 
the two-particle propagator. Numerically calculated Hubbard-I correlation function is 
presented and compared against the Monte Carlo data in Sect. [3j The results obtained 
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are discussed in Sect. HI 



2. Hubbard-I approximation 

2.1. General idea 

The idea of the Hubbard-I approximation for Hamiltonian (CD) is as follows. First, we 
introduce single-electron Matsubara Green's function G a (j — i, r) = — (TCj cr (r)c; fr (0)), 
where (...) means thermodynamic average, T is the time-ordering operator, and r is 
imaginary time. The equation of motion for G a (j\ — i, r) can be written as 

d \ 

— + /i) G a Q - i, r) = 6h6(t) + UFaa (j - i, r) 



-tVCO-iHT), (3) 



where 5y is the Kronecker symbol, <5(r) is the delta-function, and F ad .(j ~h T ) is the 
two-particle Green's function of the form 

- i,r) = -(%(r) V (r)i(0)). (4) 

Second, we write the equation of motion for function F, which includes even more 
complicated Green's functions, which, in turn, require additional equations of motion, 
etc. In the Hubbard-I approximation we avoid proliferation of these Green's functions 
by making the following decoupling: 

(fc i+aa (r)n.,(r)cl(0)) (n^)(f c j+a »cl(0)), (5) 
(rc jff (r) C J +a ,(r) Cj ,(r) C L(0)) -> (c j5 c j(T ) (6) 
x(TcJ + ^(T)i(0)) = 0, 

(rc j(T (r)cj,(r)c j+a ,(r)cL(0)) (c,J. 5 ) (7) 
x(^ i+ ^(r) C L(0)) = 0. 



To understand the nature of this approximation let us express the electron density as a 
sum of its average value and fluctuations around it: 



ru 



JO" 



(nja) + Sn^, where (8) 
Srifc = - (n Jg ). (9) 

Two other boson-like quantities, c- } -c- }a and Cj CT cj-, can be represented in the same 
manner. The above decouplings correspond to the assumption that the propagation of a 
single electron is unaffected by the fluctuations of these quantities around their average 
values. This assumption could be proven by either showing that (i) the fluctuation terms 
were small, or that (ii) there were no correlation between the single-electron motion 
and the fluctuations. Unfortunately, neither (i) nor (ii) are rigorously established. 
Because of this the Hubbard-I approximation is uncontrollable approach whose accuracy 
is unknown. 
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Using the decouplings ([MZD we derive the equation for F in the form 

d 



^(T)-f^G ff (j-i + a, r) 



This equation, together with ([3]), constitutes a closed system, sufficient for calculation 
of functions G and F. This is how the single-electron propagator is calculated within 
the framework of the Hubbard-I approximation. 

2.2. Hubbard-I approximation for the two-particle propagator 

Following Ref. [10J, we calculate here the equal-time two-electron correlation functions 
of the form 

1 N 

= E< A to A ^ + r ) + A ci(i)At(i + r )), (10) 



i=l 



where N is the number of sites on the square lattice, 

* r 

A- = ct ct 
The form-factor 

/d(r) = ^r H ,0 (5r»,l + <*r. ,-l) ~ <*r«,0 (<*r v ,l + <*r„,-l) 

corresponds to the 0^2 -wave order parameter. 

To calculate Pd(?) let us first define the propagator 

P yIm (r) = (f A i .(r)At m (0)). 

The equal-time two-electron correlation function P^(r) is related to Piji m (r) as 



-Pi,i+b,i+r,i+b'+r ( r ) 



(11) 

(12) 
(13) 

(14) 
(15) 



ibb' 



+-Pi+b,i,i+r,i+b'+r(' r ) + -Pi,i+b,i+b'+r,i+r ( T ) 
+-Pi+b,i,i+b'+r,i+r('7") 



[t -> +0, r -> — r) . 



r- »-0 



We write down equation of motion for the propagator Piji m using Hamiltonian ([I]) and 

( ^ + 2 ^) Ajlm = (-^i+aJ 1 ™ + Aj+alm) 

+5(t) ((clc la )5 im - (<w^n) + UVtsua, (16) 
where the three-particle propagator is defined as 

TW-r) = (f K(r)n Jff (r) + n i ,(r)A ij (r)j AL(0)>. (17) 
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The double occupancy is very unlikely: -Piji m = o(t/U) ~ if i = j or 1 = m. Assuming 
the absence of the magnetic order, which means, in particular, (c ia cj a ) = (q-cj-), we 
transform ( fl6|) in the Hubbard-I approximation to 

f7^: + 2Aij Pij\m=-t'^2 (-Pi+ajlm + Aj-alm) (l-^i+aj) 

^ ' a 

+S(r) (i&Mm - (c ma cl)S a ) + UV-^ m , (18) 

which is valid if i 7^ j and 1 7^ m. 

Now we need to derive an equation for V. The required calculations are onerous, 



but straightforward. They are presented in Appendix A Here we quote the final result: 



£/7\jlm« in ic ) ^(Pi+ajlm + iValm) (1-fc+aj) (19) 
a 

S{r){n ia ) ((ic lCT )5 jm - (c mCT cj a )5ii) . 

This relation is derived under the assumption that the magnetic order is absent: 
(n ia ) = (n i& ). 

Substituting ffTO]) in f[T8"j) . we obtain the equation of motion for the propagator Piji m 
in the Hubbard-I approximation valid when i 7^ j and 1 7^ m 
d \ 

^ + 2^JP ijIm = 

[-Pi+ajlm (l-^ji+a)+-Pij+alm (l-<5ij+a)] 

a 

+5(t) (1 - (n ia )) ((clcJ5 im - (c^cDSn) . (20) 

Here t = t (1 — (ni a )) if the system is uniform. 

As we wrote above, Pm m ~ o(t/U) ~ within our accuracy. Thus, jumps from site 
j = i + a to site i must be excluded. The same is true for jumps from i = j + a to j. 
To take this fact into account and to generalize the last equation for the case i = j, we 
add the term t(l - (n ia )) J2 & (^i+ajPi+ajim + <$ij+aPij+aim) to the right-hand side of ([20]). 
Therefore: 

d \ 

+ 2/i P ijlm = 



dr 

~^y^ [-Pi+ajlm (l-5ji+a-^ij)+-fij+alm ( 1 ~ ^ij+a ~ 
a 

+5(t) (1 - (n ia )) ((clcJ5 im - (c ma cl)5 n ) . (21) 

In such equation the quantity Piii m is decoupled from Piji m , i 7^ j, which is equivalent 
to the condition Pm m = 0. 

Expression (121]) corresponds to propagation of two interacting particles: the terms 
with the Kronecker symbols may be regarded as an effective repulsive interaction 
between electrons on neighboring sites with the coupling constant of the order of t(l — 
(ni a )). Thus, we reduce the strong-coupling many-body problem to the intermediate- 
coupling two-body one. 
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Figure 1. (Colour online) The two-particle propagator Pd(r), equation (fT0|) . Solid 
(red) line with right crosses is the Monte Carlo data for the doping level of 0.18 and 
U = 6 collected on 10x10 lattice, Ref. [10]. Dashed (blue) line with skew crosses is the 
Hubbard-I result for the doping level of 0.17 on 10x10 lattice. The two straight lines 
correspond to c/r 3 behaviour. 




Figure 2. (Colour online) The two-particle propagator Pd(r), equation (fT0|) . Solid 
(red) line with right crosses is the Monte Carlo data for the doping level of 0.18 and 
U = 6 collected on 10x10 lattice, Ref. [10]. Dashed (blue) line with skew crosses is 
the Hubbard-I result for the doping level of 0.18 on 28x28 lattice. Two straight lines 
correspond to c/r 3 behaviour. 



3. Numerical results 



3.1. Hubbard-I vs. MC 



In this section we present the results of our numerical calculation of the equal-time two- 
electron correlation function ( JTOj) at a given doping level. The actual calculations are 
based on Eq. ( 1C.4j) derived in Appendix B and Appendix C Doping is equal to 1 — n, 
where the average number of electrons per site is defined as n = Y2a( n ^)- The avera g e 
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Figure 3. (Colour online) Finite-size scaling for the Hubbard-I propagator: the data 
for the lattices of three different sizes (16x16, 20x20, and 28x28) and the doping levels 
close to 0.18 are plotted. While the curve for the smallest lattice (brown curve with 
right crosses) deviates substantially from the remaining two, the latter are quite close 
to each other for sufficiently small separations (log 10 r < 0.8). 



( n ia) within the Hubbard-I approximation is 

(n i(T ) = (1 - (n ia ))n , (22) 

k 

From this equation the expression for n is found: 

„ = (24) 
1 + no(fji) 

Therefore, the doping level can be adjusted by changing \i. 

For calculation at zero temperature it is convenient to use the substitution in (1C.4j) : 

^ (£) = S(^r)-" F(£) - (25) 

At T = the Fermi distribution function is zero (unity) for e > (e < 0). Contribution 
form 1/sh (s/T) vanishes at small temperature. 

Numerical results are shown in figure [Hand figure EJ where the propagator is plotted 
versus separation r = |r| on the log-log graph. The Hubbard-I propagator Pm is shown 
as the dashed (blue) line. For comparison, the same propagator computed by T. Aimi 
and M. Imada pU] using the assumption-free quantum Monte Carlo procedure, Pmc, is 
shown by solid (red) line. The subscript 'HI' stands for 'Hubbard-I', the subscript 'MC 
- for 'Monte Carlo'. 

In figure [1] we present the Hubbard-I data calculated on a 10x10 lattice for the 
doping level of 0.17. The same lattice size is used in the MC simulation of Ref. [10J. 
The MC data is collected for the doping level of 0.18. The small discrepancy in the 
densities is unavoidable in our situation. Indeed, due to the finite size of the system, the 
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number of particles cannot be changed continuously. In addition, in case of Hubbard-I 
scheme, the number of physical electrons is not an independent quantity, but rather 
has to be calculated according to formula f T2~^l) . Given these two circumstances, it is 
impossible in general to match the MC and Hubbard-I densities exactly. However, such 
a small discrepancy (less than 6%) is of little importance for our goal. 

The solid straight line fits the MC propagator in the window 0.32 < log 10 r < 0.78. 
Such fit corresponds to c/r 3 decay, as discussed in Ref. [10]. Propagator P HI is also 
fited by a straight line (blue dashed straight line) in the window 0.3 < log 10 r < 0.6. 

To extract the asymptotic behaviour of Pmif) more reliably it is better to look at 
the propagator on a bigger lattice. In figure [3] it is shown how the propagator converges 
for large lattices. We observe that the curves for the lattices with M > 20 almost 
coincide for moderate r. Thus, such systems can be used to find the asymptotic. To 
that end, let us examine figure El It presents the Hubbard-I data for 28x28 lattice 
is shown together with the MC data. Analysis of figure [2] reveals that the difference 
between assumption-free MC and much more simple Hubbard-I calculations is not large, 
moreover, qualitative appearances of the curves are similar. In particular, both curves 
decay as a power-law c/r u . In [10J it is established that v pa 3, which is consistent with 
the free-electron behaviour. To determine this exponent for the Hubbard-I correlation 
function we use two-parameter least-square fitting procedure for the Hubbard-I data in 
the range 0.3 < log 10 r < 0.9. The fitting returns the following values: 

log 10 cfit = -1.2 ± 0.1, j/piT = 3.1 ± 0.2. (26) 

The value of i^fit is consistent with the MC result. Below we will always assume that 
the Hubbard-I correlation function decays with v = 3. Deviations from the free-electron 
behaviour at large r are a manifestation of the finite-size effects. 

The fact that the MC propagator decays as the free-electron propagator is a 
surprising revelation of Ref. [ID] . The ability of the Hubbard-I approximation to 
capture the 1/r 3 law is a rather simple consequence of the equation of motion, (J2TJ) . 
Indeed, in the latter equation, if r is large, the interaction terms (the terms with the 
Kronecker symbols) contribute to the s-wave channel only; as for the <i-wave propagator, 
its equation of motion coincides, up to a normalization factor, with that for free fermions. 

In figure [2] the obtained values of the propagators are "noisy" . These oscillations 
arise, because, in general, the propagators depend not only on the distance r, but on 
the direction of the vector r as well. Despite the "noise", one can easily observe that 
the Hubbard-I propagator Phi is smaller than the Monte Carlo propagator Pmc- 

To examine Phi and Pmc m a more rigorous manner let us quantify the 1/r 3 
behaviour of the propagators. For this aim we fit both functions in figure [2] by straight 
lines in the window 0.32 < log 10 r < 0.78. The fit lines are defined by the equations: 

log 10 P MC = -0.81-31og 10 r, (27) 

log 10 P jff/ = -1.29-3bg 10 r > (28) 

which correspond to Pd{f) ~ c/r 3 with cmc ~ 0.15, Chi ~ 0.051, and the ratio 
cmc / Chi ~ 2.9. Thus, the Hubbard-I propagator's asymptotic coincides with the 
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Figure 4. (Colour online) Function [a — I(r)], equation (|29|) : solid (red) line with 
right crosses is the Monte Carlo data. Three other lines are the Hubbard-I results for 
different lattice sizes. The parameters are the same as in figure [21 Two straight lines 
fit the propagators at moderate r according to equation The values of a's are as 
follows: a M c = 4.1, ajfj = 1.45, « a™ z = 1.40. 



asymptotic of the propagator calculated with the help of the assumption-free Monte 
Carlo method within an order of magnitude, but, of course, not exactly. 

We can further compare the Hubbard-I and Monte Carlo data by eliminating 
"noise" in the propagators with the help of the following procedure. Define the sum 

J(r)= p d(*,y)- (29) 

yj x 2 +y' 2 <r 

If we assume that Pd(x,y) ~ c(<ft)r~ 3 , where <p is the polar angle, then J(r) can be 
approximated as 

r b 

I(r) ps / dxdyP(x,y) , (30) 

where 

(3=1 c(0)rf0 = 2nc, (31) 
Jo 

and the constant a can be estimated as a ps J(r max ), with r max = M/2. Function 
/(r), by its definition, is independent of vector r's direction. Consequently, it is smooth 
function, free from irregular oscillations present in figure El This makes the fitting 
procedure more robust. 

Dependence log 10 (a — J) vs. log 10 r is plotted in figure H] for the same values of 
parameters as the data in figure [2j Similarly to the data presented in figure [31 the 
16x16 curve lies somewhat away from the curves corresponding to the larger lattices. 
The curves for 20x20 and 28x28 lattices coincide with each other almost everywhere, 
suggesting that for systems of such sizes the finite-size effects are important for large r 
only. 
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Figure 5. (Colour online) Dependence of d-wave and s-wave propagators on r; Pd(r) 
- dashed (green) line and P s {r) - solid (red) line. Both propagators are calculated 
with the help of Hubbard-I approximation for the doping level of 0.19 on 16x16 lattice. 



To determine the asymptotic parameters the curves in figure H] are fitted by the 
straight lines: 

logio(« - -0 = logio 2vrc - log 10 r. (32) 

The values of the fitting constants are: cue ~ 0.15, Chi ~ 0.055. As expected, they are 
close to cmc an d c#j. 



3.2. Application 

The numerical procedure in the case of the Hubbard-I approximation is much 
simpler than the quantum Monte Carlo computations. One can use it to derive 
different properties of the Hamiltonian (J2J). As an example, we apply the Hubbard-I 
approximation (IB. 101) to calculation of the equal-time two-electron correlation function 
in the s-channel: 

1 N 

P ^ = w E< A !(i)A s (i + r ) + A S (i)Al(i + r)>, (33) 

i=l 

where 

AsW = 7f E (Aii+r + A i+ri ) , (34) 

* r 

and f s is the s-wave form-factor 

/ s (r) = 5 ry0 (5 rx i + <J r „_i) + 6 rx0 (S ryl + 8^) . (35) 

The result is shown in figure [5] by solid (red) line. Function P s (r) decays quickly for 
small r. At larger r it oscillates around finite- value plateau. As could be expected, the 
s-wave superconducting correlations are weaker and decay faster than the <i-wave ones 
[compare solid (red) and dashed (green) lines in figure [5] . This is because the electrons 
in s channel experience strong repulsion suppressing the superconducting correlations. 
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Note also that the Hubbard-I correlation functions show no sign of the 
superconducting instability in both s and d channels: both correlation functions decay 
at large r, while in a superconducting phase the order parameter correlation function 
must saturate at large distances. 

4. Discussion 

We calculated two-electron equal-time ci-wave propagator using the Hubbard-I 
approximation and compared the results with quantum Monte Carlo computations [10J. 
The results of both approaches are in good qualitative agreement. Specifically, both 
methods predict that the correlation function decays as 1/r 3 at r exceeding several 
lattice constants. This allows us to assess the reliability of the Hubbard-I approximation. 

Technically speaking, the 1/r 3 decay of the Hubbard-I propagator is a consequence 
of the fact that the interaction in the equation of motion, (12ip . is important in the 
s-channel only, whereas in the <i-channel its effect decreases when r increases. 

The Hubbard-I approximation underestimates the value of the propagator residue 
by a factor of order three. Such discrepancy is rather reasonable for an uncontrollable 
approach. 

Our results for the two-particle correlation functions suggest that the Fermi liquid 
is a poor approximation for the Hubbard model near half-filling. We show that in <i-wave 
channel the correlation function decays as r~ 3 , which is consistent with "free fermions". 
However, fast decay of the s-wave propagator (see figure [5]) does not agree with the Fermi 
liquid picture. One, however, must remember that, unlike the o?-wave propagator, the 
Hubbard-I s-wave correlation function is not compared against a controllable method, 
thus, its accuracy is unknown. 

Interpreting our calculation one should keep in mind the following restriction. Since 
the MC data for large r is not available, we cannot say how Hubbard-I approximation 
fairs for r exceeding 5-6 lattice constants. The MC data of Ref. pi)] cannot be used 
to rule out the possibility that the Hubbard model exhibits superconductivity with the 
correlation length larger than 5-6 lattice constants. Our method in its present form is 
too crude to capture the superconducting correlations. 

To conclude, generalizing the Hubbard-I approximation we derive the equation of 
motion for the two-particle propagator. For the square lattice of finite size this equation 
is solved numerically. It is demonstrated that the Hubbard-I propagator is consistent 
with the MC data up to a numerical coefficient of order 3. Thus, Hubbard-I method 
can be used to approximate the two-particle propagator of the Hubbard model, at least 
at a qualitative level at moderate separations. 
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Appendix A. Derivation of Eq. (119ft 



Here we obtain Eq. ( jl9l assuming that i ^ j and 1 / m. The equation of motion for 
^ijim is 



Aji m = (TC(r)AL(O) ) i Sir) 
where we introduced the notation: 



C 



H, (n ia + n\ a ) A. 



"jo / — Ijj 

= [(#kin + #int) , (riia + n ia )A i3 \ - 2/j,(n ia + n j(J )Ay. 
Neglecting "double occupancy" terms in the commutator with f/kin> we get 
[H kin , {n id + n j(J ) Ay] « -t ^ n i(f cf +aCT cj- (1 - 5 ji+a ) 



(A.l) 
(A.2) 



t „t 



(A.3) 



The factors of the type 1 — 5y appear for the following reason. Consider: 



n jo" C L C j+aa- 



(A.4) 



The last term corresponds to the case of two electrons at the same site i = j + a. It is 
small and can be omitted. Applying the Hubbard-I decoupling and taking into account 
that in the absence of the magnetization (ci CT C;-) = 0, we obtain 

[#kin, (n i& + n jCT )A y ] « -t^2(n^)c\ +gLa cl i 1 ~ W) 



+ (^j CT )cL C j+a ( x( 1 - 5 j+ai)- 

Then, we calculate the commutator with H int : 

[H^t, (nia + nj^Ay] = (n is + n- Ja ) [H int , Ay] 
= U(n ia + n j(7 ) 2 Ay. 

We note that n? a = n- }a and Un^rt^ Ay = 0(t/U). As a result, it holds: 

[i/int, (n iff + n j<T )Ay] « £7(n is + n j(7 )Ay. 
Collecting all terms together, we derive 

[H, (n is + n jCT )Ay] « -t(n ff ) ^ cj+a<x c t ( X ~ 5 ji+ a ) 

a 

+ C L C j+a<x C 1 - <W) + U(n i9 + ?v)Ay, 

where we denote (n^) = (n- }d ) = (n a ). 



(A.5) 
(A.6) 

(A.7) 



(A.8) 
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n is +n ia ) A ij; Aj m 
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(A.9) 



Deriving the latter formula one has to keep in mind that 5j m <5ji = 5j m 5i m = 0, when 
i 7^ j and 1 ^ m. From ( 1A.9I) it is possible to deduce: 

([(n iff +n jo ) A.., Aj m ]} « (n jCT ) (cU> jm (A.10) 

- (n») (<wc} s )5ii, 

Obtaining this result we took into account that in the absence of the magnetic order 

(Cn*<L,) = K 4) = 0. (All) 



Further, in (lA.lOj) the terms 

n i A (T c la 8 im = 0{t/U), = 0(t/£7), (A.12) 

are dropped. Besides, we performed Hubbard-I decoupling in three-site operators: 

VU CT ~ (n ia )clc la (l - «JjO, (A.13) 

c mff cLc lff cL « ((wD^j = 0. (A.14) 

In (|A.13[) we used the fact that n- ja ci a 8- 3 \ = 0. As a result, within the scope of our 
approximation, equation ( IA.1I) can be written as 

A)lm~— t( n o)/j [-Pi+ajlm (1 ~ <5ji+a) + Pij+alm (1 — 5j+ai)] 
a 

+ ([/ - 2/i) P ijlm + 5(r)(n a ) (<40*jm - (c^c^n) • 
Since P = 0(t/U), while P = 0(1), we omit V and //P terms, arriving at equation ( fT9l) . 

Appendix B. Fourier transformation 

To solve (l2~Tj) . analytically or numerically, it is useful to subject this equation to the 
Fourier transformation. Let us consider a square 2D lattice with the number of sites 
N = M x M. In Fourier space equation (j2ip reads: 



Pi 



kik 2 k 3 k4aj 



7CJ 



(e lkia +e ik2 j+2/i 



(B.l) 



Here 



aij 



+ (#ij+a + 5jj) Pij+ak 3 k 4 u;] , 



K 



1 - (nig) 
A 2 



(( C kl a C k 1( x) - ( C k 2 a C k 2 <x)) <W 3 4 2 k 4 . 
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After straightforward algebra we derive 

zt 
~~N 



S=T7 ^2 (^ki-K/ka-K/k! -q+Tks+q) Pk! -q,k 2 +q,k 3 ,k 4 ,w , (B .2) 



q 

where the number of the nearest neighbors on the lattice z = 4, and we introduce the 
notation 

7k = - ^exp(ika). (B.3) 

a 

The term 1Z in the right-hand side of equation (IB.lj) is transformed with the help of the 
formula for the single-electron Green's function in the Hubbard-I approximation: 

(c i(T c} (7 ) = -(l-(^))Go(i-j,+0), (B.4) 
where Go is the Green's function for free fermions. From this equation it follows that: 

(<wL> = (1 - - n F (e k - /i)], (B.5) 

where uf is Fermi distribution function, and is the energy spectrum. In the considered 
case of the square lattice and strong electron-electron repulsion (U ^> zt) we can write 
in the tight-binding approximation: = ztj^. 

Within the Hubbard-I approach expectation value (cj a c ia ) satisfies: 

(iO = -Mj + (l-(*. (B.6) 

This expression is incompatible with the usual Fermi anticommutation rules. This is an 
artifact of the Hubbard-I scheme, a discrepancy which is a consequence of the inexact 
nature of the Hubbard-I Green's functions. This means that the Hubbard-I Green's 
function fails at energy higher than t, or, equivalently, at |i — j| smaller than at least 
several lattice constants. 

Using the above relations we derive after simple algebra 
(l_/ n \)2 

71 = — ^ F ^ kl ~~ ^ +nF ^ k2 5klk3 5k2k4 ' ( B ' 7 ) 

Combining the expressions for S and TZ, we can rewrite ( TB.lj) in the form 

[iu + zt (7 kl + 7k 2 ) + 2/i] P kl k2k 3 k 4 u 

= ^y^(7k 1 +7k 2 +7ki-q + 7k 2 -q)-Pki-q,k2+q,k 3 ,k4,a; + 

(i - M ) 

M 2 

The propagator -Pkik 2 k 3 k 4 w is non-zero only if ki + k 2 = k 3 + k 4 . This is a consequence 
of the momentum conservation law. Therefore, it is convenient to introduce the total 
momentum k 4 + k 2 = Q and define 

Pkikfu, = • P k j) Q-k i) k /) Q-k_ f ,w- (B.9) 



q 

2 

[n F (£ki-^)+n F (e k2 -^)-l] 5 k ik 3 ^k 2 k 4 - (B.8) 



Equation flB~8|) can be rewritten as 

^k%. = E h ?^ P ^ + (B.10) 



q 
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where 



^Sq = - (7k, + 7Q-0 + 2//] 5 qki 



ft 



zt 

+ (7k, + 7Q-k, + 7q + 7Q- q ) > (B- 11 ) 

Q _(l"W) a 



" k * k / ~ M 2 

x h F (e : k l -A i )+^F(e r Q-k I -Ai)-l]4 I k r (B.12) 
Equation (IB.lOp may be used to calculate the propagator P. 

Appendix C. Computation of the equal-time correlation function 

Equation iTCHf)]) may be cast in the matrix form iu~P = hP + ft. Matrix h is symmetric 
and, consequently, can be diagonalized. We write down h symbolically as h = UdU + , 
where d is a diagonal matrix with the diagonal elements 

d = diag(^), (C.l) 

U is a unitary matrix, U + is its Hermitian-conjugated matrix. Thus, the solution of 
the problem can be expressed in the form 

P(w) = [U(iw - d)U + ] _1 ft. (C.2) 

We diagonalize numerically the matrix h and calculate the equal-time propagator 
-Pyim(+0) = (AyA} m ) performing frequency summation 

P( r ^+0) = TVexp(i™)Py (C.3) 



according the rule 

T exp (iru) (iu — d) _1 = i\T B (d), 

where N-q{e) = [e e l T — is Bose distribution function, and iV B (d) is a diagonal 
matrix: iV B (d) = diag[iV B [di)\ As a result, we derive 

P(r 0) = UA/" B (d)U + ft. (C.4) 

This equation can be used to calculate the correlation function numerically. 
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